#Forelesning 4, fredag 20.01.2017

#7# Enkel lineær regresjon, prediksjons og konfidensbånd.
mammals <- read.table("http://www.math.ntnu.no/~jarlet/statmod/mammals.dat",header=T)
attach(mammals)

logbody <- log10(body) #log(body,base=10)
logbrain <- log10(brain)
plot(logbody,logbrain)		#Lineær sammenheng?
linregr <- lm(logbrain~logbody)
summary(linregr)
abline(linregr)

#Sjekke residualene - hvordan skal de se ut?
residuals(linregr)
hist(residuals(linregr),freq=F) #Normalfordelte?
curve(dnorm(x,sd=0.30),add=T) #Hvor finner jeg sd=0.30 fra? Se summary

qqnorm(residuals(linregr)) #QQ-plott for å se om residualene er normalfordelte

plot(logbody,residuals(linregr)) #Økende residualer med økende x=logbody verdier? Eller konstant varians?
abline(0,0,lty=2)

# Estimat av sigma beregnet "for hånd"
sqrt(sum(residuals(linregr)^2)/(62-2))

#Og R^2
cor(logbrain,fitted.values(linregr))^2

# "nye" mulige x-verdier - mellom -2 og 4
predframe <- data.frame(logbody=seq(-2,4,length=100))
# konfidensintervall for regresjonslinja
confint <- predict(linregr,newdata=predframe,interval="c")	#?predict. Inspiser confint
# legg dette til som linje i plottet
plot(logbody,logbrain)		#Bedre?
abline(linregr)
lines(predframe$logbody,confint[,3],col="red")
lines(predframe$logbody,confint[,2],col="red")

# tilsvarende prediksjonsbånd - breiere eller smalere enn confint?
predint <- predict(linregr,newdata=predframe,interval="p")
lines(predframe$logbody,predint[,3],lty=2,col="blue")
lines(predframe$logbody,predint[,2],lty=2,col="blue")


#8# Illustrasjon av oppførselen til den estimerte regresjonsmodellen og  
## konfidensintervall for en framtidig forventet respons - ved simulering
## forventer ikke at dere skal skjønne alle detaljene i programmeringen nå
n <- 0
hits <- 0 #antall ganger konfidensintervallet inneholder den sanne forventningsverdien
repeat
{
  x <- 1:10
  y <- rnorm(n = 10, mean = 5 + 0.25*x, sd = 1) #Modell y=5+0.25x med sigma=1
  plot(x,y,col="red",ylim=c(2,12),xlim=c(-2,13))
  abline(5,.25,col="darkblue") # sann regresjonslinje
  
  estregr <- lm(y~x)
  abline(estregr,col="red") # estimert regresjonslinje basert på simulerte data
  xnew <- 11; ynew <- 5+.25*xnew
  ci <- predict(estregr,newdata=data.frame(x=xnew),int="c")[2:3]
  lines(rep(xnew,2),ci,col="red")
  n <- n + 1
  if (ci[1]<ynew & ci[2]>ynew)
    hits <- hits + 1
  cat("n = ",n,", hits = ",hits,", hits/n=",hits/n,"\n")
  Sys.sleep(.5) # en liten kunstpause
}

## Tilsvarende for prediksjonsintervallet
n <- 0
hits <- 0
repeat
{
  x <- 1:11
  y <- rnorm(n=11,mean=5+0.25*x,sd=1)
  n <- n + 1
  plot(x,y,col=c(rep("red",10),"blue"),ylim=c(2,12),xlim=c(-2,13))
  abline(5,.25,col="darkblue") # sann regresjonslinje
  estregr <- lm(y[-11]~x[-11])
  abline(estregr,col="red") # estimert regresjonslinje basert på simulerte data
  xnew <- 11; ynew <- 5+0.25*xnew
  points(xnew,y[11],col="blue")
  ci <- predict(estregr,newdata=data.frame(x=xnew),int="p")[2:3]
  lines(rep(xnew,2),ci,col="red")
  if (ci[1]<ynew & ci[2]>ynew)
    hits <- hits + 1
  cat("n = ",n,", hits = ",hits,", hits/n=",hits/n,"\n")
  Sys.sleep(.5) # en liten kunstpause
}


detach(mammals) # ferdig med denne data.frame'en



#1# Håndtering av kjente fordelinger i R (Dalg. kap 3)
#Repetisjon av d...., p...., q..... og r.... funksjoner
par(mfrow=c(2,2)) #Får fire panel på en side - 2 rader og 2 kolonner
curve(dnorm(x),-3,3)
curve(pnorm(x),-3,3)
curve(qnorm(x),0,1)
hist(rnorm(100))
 
curve(dbinom(x,size=10,prob=.3),0,10,type="h")
curve(pbinom(x,size=10,prob=.3),0,10,type="s")
curve(qbinom(x,size=10,prob=.3),0,1,type="s")
hist(rbinom(100,size=10,prob=.3),breaks=seq(-.5,10.5,by=1))


#2# Multippel regresjon - intro
install.packages("ISwR")
  library(ISwR)
data(cystfibr)
View(cystfibr)
attach(cystfibr)
multregr <- lm(pemax ~ age + bmp)
summary(multregr)

#1# To kategorier kodet ved en dummy variabel - sex
multregr <- lm(pemax ~ age + sex)
#plot3d(multregr)
summary(multregr)

#Full modell med alle forklaringsvariablene
multregr <- lm(pemax ~ age + sex+height+weight+bmp+fev1+rv+tlc+frc)
summary(multregr) #Hvorfor er ikke age lenger signifikant?

cor(cystfibr)

#Reduserrer modellen - hvor mange skal ut? Noe slik ... ?
multregr <- lm(pemax ~ age + fev1)
summary(multregr) 

#3# Kolinearitet, simulerte data
x1 <- rnorm(50)
x2 <- .8*x1+rnorm(50,sd=0.2)  #x1 og x2 korrelerte?
cor(x1,x2)
y <- 0.5*x1 + rnorm(50,sd=1)  #y avh. av x1, ihht. definert modell
kolin <- lm(y~x1+x2)
summary(kolin)
summary(lm(y~x1)) # Merk at betahat_1 endres når x2 tas ut - hvorfor?